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ABSTRACT 

We study the motion of test particles around a center of attraction represented by 
a monopole (with and without spheroidal deformation) surrounded by a ring, given as 
a superposition of Morgan & Morgan discs. We deal with two kinds of bounded orbits: 
(i) Equatorial circular orbits and (ii) general three-dimensional orbits. The first case 
provides a method to perform a linear stability analysis of these structures by studying 
the behavior of vertical and epicyclic frequencies as functions of the mass ratio, the 
size of the ring and/or the quadrupolar deformation. In the second case, we study 
the influence of these parameters in the regularity or chaoticity of motion. We find 
that there is a close connection between linear stability (or unstability) of equatorial 
circular orbits and regularity (or chaoticity) of the three-dimensional motion. 

Key words: Planets and satellites: rings, celestial mechanics, chaos. 



1 INTRODUCTION 

Light particles interacting with a central massive body 
are frequently encountered in diverse fields of physics. 
Electrons in rotating molecules and particular versions of 
the many-body problem in celestial mechanics are the 
most well-known examples, but some nuclear models fall 
into the same category. Increasing amounts of informa- 
tion about narrow planetary rings suggest that such rings 
are often associated with the so-called shepherd satellites 
iGoldreich fc Trcmainc_19 79). and may exist due to mech- 
anisms some what more complicated than the well known 
broad rings i Greenbe rg fc Brahic 1983 : ISmith et al. 198ll : 
ISmith et al~ 982: Mur ray et al. 20051 ). 

Recent progress in the subject suggest a generic mecha- 
nism, that does not depend on Kepler orbits, to explain the 
formation of rings around a rotating object which also holds 
for systems quite different from planetary rings. In rotating 
scattering systems, the generic saddle-center scenario leads 
to stable islands in phase space. Non-interacting particles 
whose initial conditions are defined in such islands will be 
trapped and form rotating rings. This result is generic and 
guarantees that the orbits supporting the ring structure are 
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rather insensitive to small perturbations and thus may play 
a role in different situat ions of the type mentioned above 
(Benet & Selisman 2000). So, although we are interested in 
studying planetary rings, there are many other systems that 
exhibit the same structure and therefore one can perform 
similar treatments by knowing the nature of the forces in- 
volved. 

Now then, based on the arguments mentioned above, 
one may conjecture that once the ring structure is formed 
around the rotating body, it remains stable, as revealed by 
the structure of its phase space. In this set of papers, we 
propose an analytical study about the dynamical aspects 
concerning the stability of ring structures and the relation 
with the existence of isolated islands in the phase space. In 
this first paper we focus in the effects related to the physical 
parameters like the mass of the central body, the geometry 
of the central body, the mass of the ring and the size of 
the ring. In next papers, we plan to examine the influence 
of the angular momentum of the system in its stability, as 
well as, large perturbations like the interaction with external 
satellites and their responses to resonance. 

Planetary rings consist in thin discs of cosmic dust 
and other small colliding particles revolving around a cen- 
tral planet in a flat disc-shaped region. The most spectac- 
ular example of r i ng structures are those around Saturn 
llPorco et al. 20051 : ICuzzi et al. 2010f ). but they are a com- 
mon feature of the other three gas giants of the solar sys- 
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tern; Jupiter, Uranus and Neptune possess ring systems of 
their own. Recent reports have suggested that the Satur- 
nian moon Rhea may have its own tenuous ring system, 
which w ould make it the o nly moon known to possess a ring 
system (| Jones et al. 20081 ). 

There are many possible mechanisms to explain the ex- 
istence of planetary rings, but essentially three of them are 
the most relevant: from material of the protoplanetary disc 
that was within the Roche limit of the planet and thus could 
not coalesce to form moons; from the debris of a moon that 
was disrupted by a large impact; or from the debris of a 
moon that was disrupted by tidal stresses when it passed 
within the planet's Roche limit. This last one allows us to 
predict that Phobos, a moon of Mars, will break up and 
form into a p lanetary ring in about 50 million years due to 
its low orbit (jHolsapple 200lf ). 

Additionally, we should take into account the increasing 
amount of data we get from extrasolar systems. The discov- 
ery of extrasolar planets (which up to date are more than 
300) by radial velocity measurements has provided the first 
dynamical characteristics of planets: orbital elements and 
mass. The next step will be to investigate physical character- 
istics: albedo, temperature, radius, etc, and their surround- 
ings. Among the latter are planetary rings. The emitted 
thermal infrared light from the planet should show no phase 
effect assuming the planet is in thermal equilibrium. But the 
reflected visible light will vary with phase angle, as should be 
shown by a broad-band photometric follow-up of the planet 
during its orbital motion. In particular, it has been stud- 
ied from different perspectives how the presence of a ring 
around a planet would influence its brightness as a function 
of its orbital position and based on that there have been 
multiple theoretical predictions for photometric and spec- 
troscopic signatures of rings ar ound tran siting extrasolar 
planets (Barnes & Fortncy 2004; lArnold fc Schneider 200S; 
lOhta. Taruva and Suto 20091 ). So, understanding the dy- 
namics of composed planetary systems would not con- 
stitute just an astrophysical curiosity, it would have 
some impact on the stra tegy for their detection (see 
lArnold and Schneider 20041 and the references therein) . 

On the other hand, it has been shown that collisions 
play an importan t role in the dy namical aspects of the rings 
l|Longaretti 19921 : ISicardv 20061 ). Under some physical as- 
sumptions, one can see that as this process dissipates me- 
chanical energy, while conserving angular momentum of the 
ensemble, tend to flatten the disc perpendicular to its total 
angular momentum, and to circularize the particles orbits. 
If the planet is not spherically symmetric, as it is the case 
for all the oblate giant planets, apart of being flat the conflg- 
uration becomes a equatorial ring. Of course, some physics 
is still missing in the description of the rings and an inter- 
esting question here would be if these effects that one does 
not take into account as a flrst approximation are relevant in 
the dynamics of the ring. A classical approach is to consider 
what happens to the system when small perturbations are 
applied. 

Commonly, disc-like systems such as planetary rings 
are described by a set of coupled equations: two hy- 
drodynamic equations (Euler's and continu ity equation) , 
Poisson's equation and a equation of state (|Toomre 19641 : 
iBinnev fc Tremaine 200^ ). For small enough disturbances, 
these equations can be linearized and solved under fur- 



ther simpliflcations as there are still missing analytical self- 
consistent models of planetary rings. Nevertheless, it is pos- 
sible to obtain useful information of this by using simplifled 
models which illustrates quantitatively how pressure and ro- 
tation tend to stabilize the disc against self-gravity. 

Although most rings were thought to be unstable and 
to dissipate over the course of tens or hundreds of mil- 
lions of years, recent evidence coming from NASA's Cassini- 
Huygens Mission suggest that Saturn's rings might be quite 
old, dating to the early days of the Solar system. Then, the 
exact mechanism to explain which physical factors account 
for the stability of these systems is still an open question to 
be answered. 

Our approach to study the stability of rings will be quite 
different of the mentioned above. The method we will use 
is insensible to the hydrodynamical aspects of the rings but 
with the advantage that we will have under control the geo- 
metrical aspects of the models, namely the size and shape of 
the ring, the mass quotient between the planet and the ring 
and deviations from the spherical geometry of the planet. 
Also, similar techniques can be applied to study the effects 
caused by the angular momentum of the planet and even the 
influence of external satellites, but these topics are out of the 
aim of this paper and will be left for future works. Another 
beneflt of this approach is that it becomes easy to elucidate 
a non-trivial connection between the stability in equatorial 
circular orbits and the regularity in three-dimensional mo- 
tion. In particular, we will see how the existence of isolated 
islands will lead to stable ring conflgurations and how rele- 
vant are the physical parameters mentioned above. 

Now then, the exact gravitational potential of a ring 
of zero thickness and constant linear density is given in its 
exact form by an elliptic integral that is seldom used for 
practical purposes. This potential is usually approximated 
by truncated series of spherical harmonics, i.e., a multipolar 
expansion (see for example the d eep analysis of orbits per- 
formed in iTresaco fc Ferrer 2O10l ). As far as we know there 
are not simple expressions for the exact gravitational poten- 
tial of a flnite flat ring with inner and outer edges and any 
surface density. 

Meanwhile, simple potential-density pairs for thin 
discs are known. The simpler is the Plummer-Kuzmin 
disc iKuzmin 19561 ) that represents a simple model with 
a concentration of mass in its center and density that 
decays as 1/r^ on the plane of the disc. This structure 
has no boundary even though for practical purposes 
one can put a cutoff radius wherein the main part of 
the mass is inside, say 98% of the mass. There are 
many other mo dels in the U tera ture jLvnden-Bell 19621 : 
Mestel 196^: iToomre 196^: iHunter fc Too mre 1969: 
Kalnais 1972: iJiang 200d : Ijiang fcl Moss 20oi; 

Gonzalez & Reina 2006"; Ijiang fc Ossipkov 200?! : 

Pedraza, Ramos-Caro fc Gonzalez 200^ some of them 
of infinite extension but others with an outer edge. Among 
the last ones, a family of simple models are t he Mor- 
gan and Morgan di scs |Morgan fc Morgan 19691 ). which 
has been inverted (|Lemos fc Le telier 19941) in order to 
obtain infinite discs with a central hole of the same 
radius of the original disc. We can also put a cutoff 
in these structures and therefore the inverted Morgan 
and Morgan discs can be considered as representing a 
flat rings. These models have been used by the authors 
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in order to study the superposition of an annular disk 
with a central black hole in th e context of the General 
Theory of Relativity (see also Scmcra k fc Sukova 201ol : 
iLora-Claviio. Ospina-Henao &: Pedraza 2010l . for recent 
related works). 

Another approach to construct flat ring st ructures is by 
the superposition of different kind of discs. In lLetelier" 20071 
the author superposed Morgan and Morgan discs 



Flat Rings 



l|Morgan fc Morgan 196% w hile in IVogt fc Letelier 20091 
the Kuzmin-Toomre discs (jKuzmin 19561 : iToomre 1963l l 
were used. In this paper we will use the first of the above 
models which are relatively simple and superposing the 
potential of the planet will allow us to study the effects of 
those parameters we have been talking about. 




R/a 



2 MODELS OF RINGS AROUND A 

SPHERICALLY SYMMETRIC OBJECT 

In this section we shall focus on simple analytical models 
that represent sources conformed by an axisymmetric thin 
ring and a central object that we will assume with spherical 
symmetry, i.e. a body characterized only by its monopolar 
moment. The effects concerning with oblate or prolate de- 
formation of the center of attraction, i.e. the inclusion of a 
quadrupolar moment, will be studie d in the last section. 

In a recent work, iLeteUer 20o"7l . shows how to construct 
a simple family of potential-density pairs for flat rings by 
means of the superposit ion of Morgan and Morgan discs 
l|Morgan fc Morgan 196Sj ). These discs are flnite in exten- 
sion and have a well-behaved surface mass density with a 
maximum in the center and monotonically decreasing up to 
the edge. The structures obtained by the superposition of 
discs with different densities have a finite outer radius and 
zero density on their centers, i.e. discs with a hole in their 
centers, or in other words fiat rings. Although the models 
do not have an inner edge, for practical purposes one can 
put a cutoff radius and neglect the residual density (which 
becomes smaller for higher members of the family). 

The Morgan & Morgan discs are obtained by solving 
the Laplace equation in the natural coordinates to repre- 
sent the gravitational potential of a disc-like structure, i.e. 
oblate coordinates (^, rj, if) (defined in the ranges < ^ < oo 
and — 1 < < 1) that are related to the usual cylindrical 
coordinates (R, z, tp) by 



(1) 



where a is a positive constant defining the disc radius. The 
inverse relations are given by 



ai = Re [V-R^ + - 

\sjR^ + (z - iaf 



-Im 



(2) 



Note that, according to ([T]), on the equatorial plane z = 
one has to distinguish between two regions: (i) the points 
inside the disc, with coordinates C = and r\ < ^/l — /a'^; 
(ii) the points outside the disc, with coordinates 77 = and 
C > \/i?Va2 - 1. 

The mass surface density of each disc (labeled with the 



Figure 1. The mass surface density for the first ten members 
of the ring family constructed as a superposition of Morgan and 
Morgan discs. 



positive integer n) is given by 



,2\ n-1/2 



(3) 



where M„ is the total mass and a the disc radius. Such 
mass distribution generates an axially symmetric gravita- 
tional potential, that can be written as 



"i>n(C,77) = - ^A2k,2-aq2k(S,)P2k{ri)- 



(4) 



Here, P2fc(T?) and q2k{£,) = *^*^^Q2fc(iC) are the usual Leg- 
endre polynomials and the Legendre functions of the second 
kind respectively, and A2k,2n are constants given by 



A2 



M„G7ri/2(4fc + l)(2n + l)! 



a22"+i(2fc + l)(n - fc)!r(n + k + 3/2)g2fe+i(0) ' 

(5) 

where G is the gravitational constant. If we consider discs 
of the same radius a and decreasing mass, 



M„ = 



27rEc(i 
2n + l 



(6) 



where Ec is a constant that will be taken equal for all discs of 
the Morgan and Morgan family, we obtain discs with surface 
density 



2\ 1/2 



1-^ 



2 \ 



(7) 



In order to obtai n a mass dens ity in agreement with a flat 
ring distribution. iLetelier 200?! considered the superposition 

\n — k 



k=0 



Ec 1 



R^\ R^" 



(8) 



We have that the above superpositions give rings of radius a 
and a central residual density that becomes smaller for larger 
n. For practical purposes one can put a cutoff inner radius 
bn, which we will take as the radius such that the density 
falls below the 1% of its maximum. In table [1] we show the 
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n 


bn/a 


n 


b„/a 


1 


0.06210 


6 


0.59222 


2 


0.23292 


7 


0.64335 


3 


0.36991 


8 


0.67884 


4 


0.43940 


9 


0.70798 


5 


0.54302 


10 


0.73231 



Table 1. Values of the ratio bn/a, for the first ten members of 
the family of rings given by the mass surface ((Sj. Here b ji repre- 
sents the inner cut-off radius and a is the outer radius. For higher 
members of the family the mass concentration tends to be located 
near the outer edge. 



values of 6n/a for the first ten members of this family and, 
in figure [T] their corresponding mass surface densities. The 
total mass of each ring is 



M„ = 



27rEc 

3/2 2 



1 _ 4 ] r^"+Ur, 



1/2 



2r(n + 5/2) ■ 



(9) 



Thus by increasing n, the mass of the flat ring decreases 
along with the size of the region where it is distributed (near 
the outer radius). Note that this is a feature associated to the 
particular family of models we are considering, not with the 
corresponding physical applications. The values of n and Ec 
are inferred from basic data corresponding to a special case. 
For example, to describe a ring with mass 6 x 10^'Kg, mner 
radius 122000Km and outer radius 137000Km, we have to 
set n = 28 and Ec ~ 19>QQg / cm? , leading to to a maximum 
mass density of about AQg/crn^ (these values agree roughly 
with the measurement s of the so-called Ring A of Saturn 
jPoughertv et al. 20091 ')'). 

The potentials associated to these structures can be ob- 
tained using a superposition with the same coefficients as the 
ones used for the densities, i.e.. 



»!(-!)"-' 
{n-k)\k\ 



' rt ' 



(10) 



where $Jj is the same as "1>„ with the masses given by equa- 
tion (O. 

Finally we add a monopolar term, which represents the 
exterior field of a central spherically symmetric object, to the 
ring potential described above. Thus, the total gravitational 
potential reads 



GMr, 



(11) 



where Mp is the mass of the central object. Although our 
4-parametric toy model is quite simple, we believe that this 
is a good starting point to study the effects on the dynamics 
caused by the mass ratio M^/Mp and the radius ratio a/b^ 
which is one of the purposes of this paper. The effects caused 
by the rotation of the central body (its own angular momen- 
tum would be an additional parameter), will be studied in 
a next paper through the post-Newtonian scheme. 



3 MOTION OF TEST PARTICLES 

Let us deal with the problem of motion of test particles 
around the models described above. Since $„ is static and 
axially symmetric, the specific energy E and the specific 
axial angular momentum £ axe conserved along the particle 
motion, which is restricted to a three dimensional subspace 
of the {R, z,Vr,Vz) ph ase space. Each orbit is d etermined 
by the set of equations (|Binnev fc Tremaine 2008h 



R = Vr, z = V„ 

oR OZ 
where $^ is the effective potential, defined by 

*2 



2i?2' 



(12a) 
(12b) 

(13) 



in terms of which the particle's energy can be written as 

(14) 



Relations (|12ap - (|12b|) define an autonomous system whose 
equilibrium points are Vr = = z = Q and R — Rc, where 
Rc must satisfy the equation 



OR 



(flcO) 



i?3 



+ 



9£n 

dR 



= 0, 



(15) 



(Ra,0) 



that is the condition for a circular orbit in the plane z = 0. 
In other words, the equilibrium points of the system occur 
when the test particle describes equatorial circular orbits of 
radius Rc, specific energy E = $*(i?c,0), and specific axial 
angular momentum given by 



Rf- 



OR 



(16) 



(HcO) 



The description of circular orbits in the equatorial plane is 
a first step to understand the linear stability of the system, 
as well as, the regularity or chaoticity of three dimensional 
orbits. On the one hand, if one assumes as a first approxima- 
tion that the structures are built from particles moving only 
in circular orbits, the epicyclic and vertical frequencies of 
quasi circular orbits provide us a criterion for the system's 
stability (see section [3]). On the other hand, the analysis 
of such frequencies also leads to determine the existence of 
saddle points, which are preliminary indicators of irregular 
motion. 

In order to deal with the problem of correlation between 
regularity of three dimensional motion and the stability of 
circular orbits, we have to distinguish between exterior and 
interior motions of test particles, separately. As we pointed 
out in the last section, the relation between oblate and cylin- 
drical coordinates is given by ([2]) or, in a more explicit form 



^ = A{R,z) cos [e{R,z)], 
rj = A{R,z) sin [e(R,z)] , 



where 



A{R,z) = 



Q{R,z) 



42" 



+ 



R^ + z^ 



- 1 



1/4 



■ arctan 



2za 



R2 + . 



(17) 



(18) 



(19) 
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for points located outside the disc zone, while for points 
inside the disc region, we have 



Tj = ±^1 - 7?2/a2, for z = and < R < a. 



(20) 



The ambiguity in the sign in the last equation is due to the 
singular behavior of the coordinate rj on crossing the disc. 
Rigorously speaking we have rj = \/l — /a^ at 2: = 0^ 
(upper limit) and rj — — y^l — W- jo? at z = 0~ (lower limit). 

For later formulae, it is important to point out that 
the piecewise form of this transformation rule leads also to 
a piecewise form in the first and second derivatives, when 
evaluated at the equatorial plane. After some calculations, 
we obtain the following relations for the first derivatives of 
^ and 77, at 2 = 0: 



-1/2 



d^/dR = dv/dz = 0, 



for < 7? < a. 



dv/dR = T{R/a-){a^ - i^')"'''^ 
and 

d'n/dz = {R^ - aY^I'', d-q/dR = d£,/dz = 0, 



di/dR={R/a){R^ 



for R> a. 



(21) 



(22) 



The ambiguity in the sign in (|21|l has the same meaning 
as in (|17p. These relations are used to compute the second 
derivatives in the equatorial plane and the result is 



d^r]/dR^ = =Fa(a - R ) 
d^ri/dz^ = ±{R^/a){a^ ^ 



-3/2 



^2.-3/2 



for 0< R<a, 



(23) 



and 



= -a(7?' 



a ) 



-3/2 



d^^/dz^ = {R^/a){R^ 



-3/2 



for R> a. 



(24) 



The latter two equations will be important in the calculation 
of epicyclic and vertical frequencies (see section Q. 

By introducing (|2ip in the equations of motion (|12b|) . 
when evaluated in the inner zone, we can derive the equa- 
tions of motion for a test particle in the case in wich z — 0, 
R < a. For Vr we have 



Vr 



i?3 



GMp 
R^ 



R 



Wa'^ - R^ 9t] 



5=0 



(25) 



Note that, in the derivation of this equation, there is no 
difference between the choice — \J\ — R? jo? along with 
drj/dR — {R/a){a? — R^)~^^'^ , and the other option, 77 = 
-^/T^^WJa? along with dn/dR = -{R/a){a^ - i?^)-^/^ 
This is due to the fact that d'^n/dr] is an even function of 
r]. In contrast, the equation for 14 has a change of sign on 
crossing the disc: 



5=0 



1 9£n 

- i?2 di 



5=0 



for 2 = 0+ (26a) 



for 



0". (26b) 



However, since has symmetry of reflection with respect 



to the plane z — Q, its ^-derivative must vanish exactly in 
the equatorial plane and we have 



v. = 0, 



for 2 = 0, 



(27) 



ensuring the existence of circular orbits inside the disc. Now 
then, according to (|25|) . one can verify that equation (|16p . 
for inner equilibrium points, can be cast as 



Rt 9£n 

a^a'^ ~ R'i d-q 



«=o 



for Rc < a. 



(28) 



On the other hand, for outer circular orbits (equilibrium 
points outside the disc) equation p6|) becomes 



Rt 



9C 



?7 = 



5=yfli/a2_l 

for Rc > a. 



(29) 



Equations (|28p and (|29p are relevant in the derivation of 
the quadratic epicyclic and vertical frequencies, in order to 
provide a criteria for the linear stability of the structures 
studied here. 



4 LINEAR STABILITY OF STRUCTURES 

As it was mentioned above, we assume the simplified model 
that the structures are built from particles moving in con- 
centric circles. In this first approximation, the stability anal- 
ysis of circular orbits associated to test particles provides a 
stability crite rion for the structure , assumed to be a rotating 
ring of fluid jLord Ravleigh 19161 : iLandau fc Lifshitz 1987[ 
iLetelier 20071 ). For this reason, we now examine the behav- 
ior of the epicyclic and vertical frequencies associated to 
quasi circular orbits. These quantities describe the response 
of test particles to radial and vertical (2-direction) perturba- 
tions, when describing a circular motion. The epicycle fre- 
quency K and the vertical frequency v, can be calculated 
from the effective potential "3?^ through the following rela- 
tions (Binnev fc Tremaine 200&1 : 



dR^ 



(RcO) 



dz^ 



(30) 



If relation H16|) is introduced in the second derivatives of 
$* , we can obtain and as functions of R^. Thus, val- 
ues of Rc such that k > and (or) v > corresponds to 
stable quasi circular orbits under small radial and (or) ver- 
tical perturbations, respectively. Otherwise we find unsta- 
ble circular orbits. Since the equatorial plane is composed 
by two regions, i.e. inside and outside the disc, we must 
define each of these quantities as piecewise functions. The 
quadratic epicyclic frequency for equilibrium points in the 
inner zone can be written as 



7?3 



ml - 4a^ 



5=0 



a2 (a2 - Rl) drj^ 



for Rc < a, 
(31) 
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whereas that for equihbrium points in the outer zone we 
have 



100 



2 _ GMp ml - 4a^ 

Rl 



77=0 



+ 



for Rc > a. 



(32) 



Similar expressions can be derived for quadratic vertical fre- 
quency, which inside the disc takes the form 



GMr, 



+ 



Rr 



Rl a{a^-Rlf'^ drj 



+ 



a? - Rl de 



for Rc < a, 



(33) 



and outside the disc is 
2 GMp , 



Rl 



+ 



Rl aj'n 

a(i?2_„2)3/2 



r|=0 



Rl — drf 



for _Rc > a- 



(34) 



Equations (|31|l and (|33|l are used to determine the conditions 
to be satisfied by the parameters of the system, so that the 
stability condition is fulfilled (the relations (|32|) and (|34p are 
studied in the next section). By keeping Ec and a constants, 
we can search the set of Mp values that make stable configu- 
rations. From here on we shall use G = Ec = a = 1 without 
loss of generality. In Figure [2] at the left side, we can see 
the behavior of (inside de disc) for different values of the 
planet's mass, in the case n = 1. Since it is a positive con- 
cave function with a critical point in the range < -Rc < a, 
we can find the minimum value of Mp (or the maximum rate 
Mi/Mp) for which {Alp, Rc) is positive in this range. We 
find such value by solving the simultaneous equations 



U(M;,i?:) = o, 



k\m;,r*c)^o, 



for the variables M* and R*, representing the minimum 
value of the planet's mass and the corresponding critical ra- 
dius, respectively. We can see that by increasing Alp, starting 
from Mp, we obtain increasing values for k^. 

Figure [3] shows the behavior of the quadratic vertical 
frequency, also for the model n = 1, and we note that it 
is a monotonically decreasing function, in the range < 
-Rc < a, with a minimum at -Rc — a. In this case, to find 
the minimum value of Alp (which we shall denote as Mp*) 
for which {Alp, Rc) is positive, it is enough to solve the 
equation 



v\ai; 



,a) = 0. 



for the variable M" . As in the previous case, we note that by 
increasing Alp, starting from M** , the values of increase. 

Another feature showed in figures [2] and [3] is that, for 
large values of Mi/Mp, the behavior of and in the 




Rr 



Figure 2. Behavior of quadratic epicyclic frequency for the model 
n=l, for Mi/Mp = 0.05, 0.4, 0.8, 1.2, 1.6, 2, 2.8. For values larger 
than Mi/Mp = 2.8 (bottom curve on the left and right sides), is 
negative in the equatorial plane. For smaller values, is greater 
and the gap of discontinuity at Rc = a = 1 is smaller. 
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Figure 3. Behavior of quadratic vertical frequency for the model 
n = 1, for Nh/Mp = 10"", 0.01, 0.02, 0.0485. For values larger 
than K'I\/Mp = 0.0485 (bottom curve on the left side and top 
curve on the right side), is negative in the equatorial plane. 
For smaller values, u'^ is greater and the gap of discontinuity is 
smaller. 



range -Rc > a is very different from the behavior in < 
Rc < a. In general, we see that 

2 



lim 



2 _L 



lim 



(both are finite values), but the difference between these two 
limits is attenuated by decreasing the ratio AIi/Mp. We say 
that there is a gap of discontinuity at -Rc = a. The same 
statements hold for and models with n > 2. It is clear 
that the gap disappear as AIi/Mp — >■ 0. 

The behavior of epicyclic and vertical frequencies 
sketched above is the same for the other models n = 2, 3, . . ., 
so that we can compute the maximum values of the rate 
Mn/Alp leading to > or > 0. Such values are listed 
in Table (2] for the first ten members of the family (|lip . 
Since in all cases Mn/Mp* < Mn/Mp, we can establish that 
Mp* is the parameter that determines the boundary between 
stability and instability in configurations characterized by 
gravitational potentials of the form (|11|) . 
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n 


Mn/M; 




n 


m„/m; 


Mn/M;* 


1 


2.83102 


0.04850 


6 


0.13148 


0.00395 


2 


0.86297 


0.02141 


7 


0.10024 


0.00303 


3 


0.43398 


0.01205 


8 


0.07909 


0.00240 


4 


0.26554 


0.00772 


9 


0.06407 


0.00194 


5 


0.18068 


0.00537 


10 


0.05299 


0.00161 



Table 2. Ratios M„/M* and Mn/M;* for the first ten members 
of the family of configurations represented by Ulljl . Here M* and 
M;* are the minimum value of the central body's mass such that 
> and i/^ > 0, respectively. In all cases Mn/M; > M„/M**. 



0.05 - 



Staibilit J of the Monopole -Ring System 




Figure 4. Correlation between the rate Mn/M" and bn/a ac- 
cording to the behavior of quadratic vertical frequency and 
quadratic epicyclic frequency . The interpolation line among 
points is a separatrix between a stability grey region, where 
> and > 0, and an instability white region, where <0 
(and K? may be positive or negative). The points plotted here 
correspond to the first ten models, from left to right. 



It is important to note that, in these models with fixed 
exterior radius a, the smaller the size of the n-th ring model, 
the larger is M** . The reason is that, for higher members 
of the family (jlip . the ring's mass concentration tends to be 
located near the outer edge and, therefore, away from the 
central monopole. Thus, it will be required increasingly cen- 
tral body mass to provide stability to the ring structure, as n 
grows. This fact can be glimpsed through a comparison be- 
tween tables [2] and [T] from which one might infer that there 
some kind of correlation between Mn/M** and the ring's 
size of the family of models. In figure |4] we show this corre- 
lation by plotting the values of tables [2] and [1] and interpo- 
lating the corresponding points, for the first ten members of 
the family. Thus we can sketch a boundary between stability 
and instability of self-gravitating configurations under study. 
In figure 21 points located in the grey zone corresponds to 
parameters leading to stable configurations, while the points 
in the white region are associated to unstable configurations. 
We find that the points belonging to the sepparatrix of fig. 
[4] can be fitted by the relation 



Mn/M** ^ 0.065exp(-4.676„/a) . 



(35) 



5 THE PHASE SPACE STRUCTURE 

In this section, we shall make a description of three- 
dimensional motion through the study of phase space struc- 
ture associated to orbits of test particles. We are principally 
interested in the influence of the mass ratio, in the regu- 
larity of three-dimensional bounded orbits. The influence 
of the size of the ring can be inferred from it, because in 
the models studied here, fe„ is automatically determined by 
Mn- This influence has already been analyzed in relation 
with the stability of equatorial circular orbits and now we 
extend such study to the case of more general orbits (re- 
member that the motion restricted to the equatorial plane, 
which is completely integrable, can be exclusively classified 
as stable or unstable). In order to show how the nature of 
bounded motion is conditionated by the linear stability of 
the self-gravitating structures, we will use parameters close 
to the critical values we have defined previously (table [2]). 
The infiuence of the quadrupolar deformation is considered 
in section [S] 

We shall use cylindrical coordinates and plot 2 = 
surfaces of section corresponding to the equations of 
motion H12ap - (|12b|) . for different values of the conserved 
quantities E and £. It is worth to point out that we need 
to make explicit distinction between the orbits that cross 
the 2 = plane at < J? < a and the ones that cross it 
only a.t R > a. The former are called disc-crossing orbits 
(DCO) and the latter we shall denote as non disc-crossing 
orbits (NDCO). The reason is that, for each of the above 
situations, the origin of irregular motion is different. For 
the case of DCO the discontinuity in the z-component 
of the gravitational field (equations (|26ap - (l26bp ) can 
produce a fairly abru pt change i n the cur v ature, leading 
to irregular motion JSaa 19991: ISaa 2000l: iHunter 2003 : 
iRamos-Caro. Lopez-Suspez fc Gonzalez 2008 ). Somehow, 
this is a problem analogous t o the case of the chaotic behav - 
ior of Chua's circuit (|Matsximoto. Chua fc Komuro 19851 ). 
which is described by an autonomous system of the type 
X = /(x), where /(x) is a piecewise function of class C° 
(continuous but not differentiable). It is the first example 
where the existence of such class of function leads to the 
existence of a chaotic attractor in a dynamical system 
(|Madan 1993h . 

In contrast, in the case of other three-dimensional or- 
bits, chaotic motion is due to the existence of saddle points 
in the effective gravitational potential. Such saddle points, 
in the particular case of a potential as will be located at 
equatorial plane, outside the disc (note that it is not possible 
to deflne saddle points inside the disc, due to the disconti- 
nuity in the potential's 2-derivative) . 

Equations (I32|l and (|34p help us to investigate the ex- 
istence of saddle points in the potential $*, through the 
evaluation of the quantity A = {d^$'^/dR^){d'^^/dz^) - 
(d^^^/dRdz)^ , which is negative when evaluated at saddle 
points. Since is symmetric with respect to 2 = 0, the term 
d^^n/dRdz vanishes at equatorial plane and the condition 
for existence of saddle points is reduced to A = k^v^ < 0. In 
the right side of figures [2] and [3] we can observe the behav- 
ior of and i/^ < 0, respectively, and the product between 
them is showed in flgurefS] for the case of model m = 1. For 
different values of AIi /Mp we note that there is a region of 
saddle points closed to the outer edge, even for the maxi- 
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Figure 5. Behavior of the product k^v'^ for the values 
Mi/Mp = 10-^,0.01,0.02,0.0485 (curves from top to bottom). 
For Mi/Mp = 0.0485, the maximum value leading to a stable 
configuration, there is a small region of saddle points near to the 
outer edge (right side of the dashed line). This region becomes 
smaller as Mi/Mp decreases. 



mum value 0.0485 leading to a stable configuration. When 
the ratio Mi/Mp decreases and we have even more stable 
configurations, the range of saddle points decreases as well 
as the gap between the values of k^i/^ near the outer edge. 
This suggest that more and more stable configurations, pro- 
duce less and less irregular orbits. 

In order to illustrate the above statements, we plot the 
surfaces of section corresponding to motion around the first 
two models, by choosing several values for the parameter 
Mp. We solve the equations of motion (|12a[) - (|12b|l by the 
Runge-Kutta 4th-method with variable time step and in- 
corporating the Henon algorithm (.Henon 19 82). in the LCP 
(Laboratorio de ComputaQao Paralela) at IMECC. In the 
algorithm we take into account explicitly the discontinuity 
in the field force by implementing the set of equations (|26ap - 
H27p . as well as, the piecewise transformation H17p - (|20p . We 
choose initial conditions at z — 10^^^, Vr — and several 
values for 7? (the component 14 is given by (|14|) and does 
not vanish). Orbits were integrated with a precision char- 
acterized by a relative error of 10"* or less, in the energy 
conservation. Integrations were carried out on times of the 
order of 10^ 

In figure[H]we have set Mi/Mp = 2.8, which determine a 
radially stable configuration, but vertically unstable. There 
is a variety of central K AM curves and small resonant islands 
outside the disc zone as well as in the inner region. They 
alternate with two chaotic regions, one of them is due to 
DCO (the most prominent) and the other, more dense, is 
the result of two orbits near the saddle point {Rc = 1.014 for 
this case). Here we have a situation where an unstable model 
admits a variety of irregular orbits. If we choose a more 
stable configuration, i.e. characterized by parameters near 
the sepparatrix of figure [l] we obtain surfaces of section as 
the shown in figure[7|where the chaotic region is restricted to 
the small central annulus enclosing the three KAM curves on 
the rigth. Note the presence of two chains of resonant islands 
enclosing the stochastic zone and other two small chains 
within it. In this case we find no chaotic motion associated 
to disc crossings, but only to exterior saddle points. 




-lii - 




-"-RSI 




-1.0 — 

0.0 1.0 2.0 3.0 4.0 

R 

Figure 6. Surface of section for some orbits around the model 
n = 1, with Mi/Mp = 2.8. Such a rate determine an unstable 
situation, where ft^ > (marginally) but < 0. We find two 
chaotic regions: (i) A prominent zone due to DCO and (ii) a 
smaller zone corresponding to two orbits with initial conditions 
next to the saddle point. In this case, we have chosen E = —0.3 
and e = 0.32. 
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1.0 



0.5 



0.0 



-0.5 



-1.0 



-1.5 





R 

Figure 7. Surface of section corresponding to three-dimensional 
orbits around the model n = 1, setting Mi/Mp = 0.0578. The 
model represents an unstable situation (although very near to 
the separatrix of figure |4jl and there are some chaotic orbits. The 
irregularity associated to these motions is due to the proximity to 
saddle point near to the disc edge. The orbits plotted have values 
E = -6 and £ = 3.6. 
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Figure 8. By maintaining tlie same values for E and £ as in 
figure [7] but decreasing Mi/Mp to 0.0485, we find a situation 
when the n = 1 model is stable (also very near to the separatrix 
of figure UJl and we note only regular orbits. 



When we turn to the other side of the border, toward 
the region of stable configurations in figure [J] regularity is 
the most common feature of three dimensional orbits. Such 
is the case illustrated in figure [S] where we find only smooth 
KAM curves along the surface of section, even in the case 
of DCO orbits. In this case we can see a chain of 29 small 
resonant islands enclosing another chain of 3 central islands. 
This transition from chaos to regularity, characterized by the 
apparition of increasingly number of islands when we pass 
from unstable to stable structures, can be seen in figures [9] 
and 1 101 which correspond to the model n = 2. In the former 
there is a sequence of three prominent central islands, one 
of them enclosing a small chaotic region with a chain of 17 
small resonant islands. In figure [10] the stochastic region is 
absent and we see only a dashed KAM curve enclosing three 
large islands. 

It is clear that we expect a similar behavior for the 
remaining members of the infinite family: the increment of 
the mass of the spherical central body favors the existence 
of isolated islands in the phase space. 



6 ADDITION OF QUADRUPOLAR TERM TO 
THE CENTRAL BODY 

It is known that in the universe there exist many centers 
of attraction with a certain deviation from spherical sym- 
metry, so we are interested in including a quadrupolar term 
in the description of the central object that makes up the 
structures studied above and examine the effects of deforma- 
tion (preserving the axial symmetry) . Now the gravitational 
potential of the structure is 



(/3) 



2(i?2 + ^2 



2.0 



1.0 



0.0 



-2.0 




R 



Figure 9. Surface of section corresponding to three-dimensional 
motion around the model n = 2 with M^jM^ = 0.0340, for orbits 
with E = —6 and £ = A. Since we are dealing with a configuration 
belonging to the unstable region of figure |4l there is a stochastic 
zone between a variety of regular KAM curves. 
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2.0 



^ 0.0 
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Figure 10. By maintaining the same values for E and £ as in 
figure (9] and choosing M2/Mp = 0.0214, we obtain a situation 
when the n = 2 model is stable and we note only regular orbits. 



where j3 is the quadrupolar moment that quantifies the 
oblate {j3 < 0) or prolate (/3 > 0) deformation of 
the central body. In general, it is related to the mass 
dens ity p(r, 6) (axially sy mmetric) through the equation 



5/2 ' 



(36) 



(tBinnev fc Tremaine 2008 

I3 = 2tv r*dr d0 sin eP2{cos 6) p{r, 9), 
Jo Jo 



(37) 
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where r = \/R^ + z'^ and cos6' — z/^/R^ + are the spher- 
ical coordinates. 

The equations of motion for test particles around such 
configurations are the same p2al) - (|12b|) but replacing $„ by 
\ Then, the relation for Vr inside the ring becomes 

5=0 



100 



, 3/3 GMp 



R a-i>„ 



(38a) 

while the relation for 14, in the inner, remains the same as 
in section [3l As a consequence, the relation that determines 
£c, for inner circular orbits, changes to 



2Rc aVa^ - R'i dr] 



for Rc < a, 



and, for outer circular orbits 

£^-GM,R.- ^^^ + 



7J = 



for 7?c > a. 



Therefore, the relations for the quadratic epicyclic frequency 
turns to 



2 _ GMp 6P_ ml - 4a' a$n 
' R^ ^ Rl^ a{a^-Rlf/' dv 



5=0 



a2 (a2 - _R2) ar;2 



5=0 



>7 = yi-H2/a 



and 



2 GMp 6^ 3i?^_4a2 a$„ 

«^ = — J::;^ ^- + 



ii3 iJ5 a(^2_„2)3/2 



for Rc < a, 
(41) 

?7 = 



+ — 



Rl d^^„ 



a2 {R2 - a2) 



77=0 



5=V-Ri/''"-i 
while for the vertical frequency we have now 



5=V-Ri/''^-i 

for Rc > a, 
(42) 



2 GMp 9/3 , 



Rf. 



R3 2RI a{a^-Rlf'' dii 
5=0 



5=0 



+ 



a2 - Rl 



„=0~r|7^ 



for i?c < a, 



(43) 



and 



2 GMp 9/3 , 



Rc 



7?3 2i?5 a(7?2-a2)3/2 9? 

>7 = 



?7 = 



+ 



R1 — 9r;2 



5=7^17^ 



5=Vfli/""-i 
for _Rc > a. 



(44) 



We note that the epicyclic frequency behaves in a very 
similar fashion as in the previous case of a spherically sym- 
metric central body. This can be seen in figures [TTI and [T2l 




Rc 

Figure 11. Behavior of for M\/Mp = 0.05 and /3 = 
0.1,0.2,0.3,0.4 (curves from top to bottom). Maintaining the 
central monopole fixed, the epicyclic frequency grows as the 
quadrupolar moment decreases. 



100 




Rr 



Figure 12. Behavior of k? for /3 = 0.1 and Mi/Mp = 
0.05, 0.15, 0.25, 0.35 (curves from top to bottom). Maintaining the 
quadrupolar moment fixed, the epicyclic frequency grows as the 
central mass increases. 



In the former we have chosen a certain fixed value for Mp 
and plotted k' for some positive values of quadrupolar mo- 
ment (assuming prolate deformation), while in the latter 
we have fixed P and plotted for different values of Mp. 
These figures reveal that the quadratic epicyclic frequency is 
a positive concave function with a critical point in the range 
Q < Rc < a. Needless to clarify that it is a monotonically 
decreasing function if we use negative values for /3. Also, we 
have to point out that this quantity has the same behavior 
for the remaining models n = 2,3, ... (the same statement 
holds for the vertical frequency). 

In contrast, the addition of a quadrupolar term intro- 
duces new features in the behavior of the quadratic vertical 
frequency. It can be observed from figures [TH] and [2] that 
v' is now a negative concave function with a maximum in 
the range Q < Rc < a, for positive values of quadrupolar 
moment (we have used the same values for the parameters 
as in figures [TTI and [T2 |l . It is worth clarifying that, by us- 
ing negative values for /3, the quadratic vertical frequency 
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Figure 13. Behavior of for Mi/Mp = 0.05 and /3 = 
0.1, 0.2, 0.3, 0.4 (curves from top to bottom). Maintaining tlie cen- 
tral monopole fixed, the range where the quadratic vertical fre- 
quency is positive, grows as the quadrupolar moment decreases. 
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Figure 14. Behavior of u'^ for /3 = 0.1 and Mi/Mp = 
0.05, 0.15, 0.25, 0.35 (curves from top to bottom). Maintaining the 
quadrupolar moment fixed, the range where the quadratic vertical 
frequency is positive, grows as the central mass increases. 



becomes a monotonically decreasing function, as in the case 
of a spherically symmetric central body. 

According to the above statements, the task to find the 
limiting values of Mp and /?, leading to stable configurations, 
is reduced to formulate the simultaneous equations 

i/'(M;,/3t,6„) =0, iy\Ml,l3\a)^0, 

for the variables and p\ i.e. the minimum value of 
the central body's mass and the maximum value of the 
quadrupolar moment such that the quadratic vertical fre- 
quency (and evidently k^) is positive in the range between 
the cut-off radius and the outer radius. In table [3] we have 
listed the corresponding values of these quantities for the 
first ten members of the family. Consequently, we show the 
correlation between the logarithm of /3n//3^ and Mn/Mp, 
plotting the separatrix between the stability and instabil- 
ity (fig. lisp . Here, Pn represents the quadrupolar moment 
of the n-th ring model and it is computed using equation 
(|37|) . We find that the variables Xn = logj^Q(|/3„|//3) and 



n 


/3n//3t 


m„/a4 


n 






1 


-16.1076 


0.04804 


6 


-0.01300 


0.00253 


2 


-0.55925 


0.02024 


7 


-0.00807 


0.00177 


3 


-0.12408 


0.01040 


8 


-0.00537 


0.00129 


4 


-0.04713 


0.00602 


9 


-0.00375 


0.00097 


5 


-0.02301 


0.00378 


10 


-0.00273 


0.00074 



Table 3. Ratios Mn/Mp and /3„//3t for the first ten members of 
the family given by II36I I. Here Mp and /3t represents the minimum 
value of mass and maximum value of quadrupolar moment such 
that the structure is linearly stable. 

Cuadrupokr Moment -Mass CorrektiDn 




-1.0 -0.5 0.0 

Figure 15. Correlation between the rate M„/Mp and Pn/li- The 
grey region represents the set of values leading to stable struc- 
tures, while the white region corresponds to the stable ones. 

Un = log]^Q(|A/„|/A/j5) can be fitted by the relation 



-0.099x„ -f 0.3452;„ - 1.592, 



(45) 



providing an approximate expression for the separatrix of 
figure [T5] 

As it was shown in lGueron fc Letelier 200 ll . for the case 
of monopole-quadrupole configurations, the chaoticity is in- 
duced by prolate (rather than oblate) deformation. In the 
case studied here, we have also a ring structure with oblate 
shape and one would expect two situations: (i) an "attenu- 
ation" in the chaoticity, for the case of a central body with 
prolate deformation; (ii) only regular motions, if the central 
body has also oblate shape. 

The surfaces of section corresponding to structures lo- 
cated at stability region of figure [TJ] are not very different 
from those shown in plotslSlor llOl However, some differences 
appear when we deal with unstable configurations. For ex- 
ample, the outer stochastic region due to DCO, in figure [S] 
disappears when we turn on the quadrupolar moment, but 
the inner chaotic region is now most prominent and shifted 
slightly to the left (details are shown in figures [TBI and [T71) . 



7 CONCLUDING REMARKS 

Simple models as the ones described by relations ^ and 
(|10l) . provide us an useful tool to analyze and understand 
the dynamics of astrophysical objects that can be modeled 
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Figure 16. Detail of the central part of figure|6]corresponding to 
n = 1, Mi/Mp = 2.8E = -0.3 and £ = 0.32. The inner stochastic 
region is generated by two orbits with initial conditions near the 
saddle point. The outer stochastic region is due to a DCO. 
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Figure 17. We have used the same parameters and initial con- 
ditions as in figure [T6l but switching on the quadrupolar moment 
to /3 = 0.03. The stochastic regions associated with this unsta- 
ble configuration are distorted, with respect to the figure [T6l The 
outer chaotic zone of DCO is absent, while the inner chaotic zone 
is more prominent. 



as a spherical (or quasi-spherical) mass surrounded by a flat 
ring structure. One important and fundamental aspect in 
the dynamics is the stability against proper modes, if we 
consider in a first approximation that the ring structure is 
formed principally by particles moving in c oncentric circu - 
lar orbits. It is known, in previous studies ()Letelier 2007l l. 
that the structure of one ring standing alone has limited 
stability but, by adding a central monopole, it can improve. 
In this study we were able to verify such statement and 
find the set of values for the parameters leading to linearly 
stable configurations. Interestingly, we find monopole-ring 
configurations belonging to the stability region of figure HI 
when dealing with parameters of the order of physical mea- 
surements performed in the solar system. For example, a 
structure with the dimensions of Saturn has Rmin/Rmax ~ 
0.1 (ratio inner-outer radius) and Mring / M^ianet ~ 10~*. 
Other example is a structure with Rmin/Rmax ~ 0.4 
and Mri na / M^ianet ^ lO"^'^, similar to the par ameters of 
Jupiter (|EUis et al. 2007I . iDoughertv et al. 20091 '). We have 
to point out that a deeper stability analysis must include 
the effect of mic rostructure in the ac t ual planetar y rings 
(see for example iMerlo fc Benet 20071: [Murray ct al. 20081 : 



iBenet fc Merlo 20091 : ICharnoz 20091 : ICharnoz et al. 2010l) 

On the other hand, there is a close connection between 
the stability of the configurations studied here, and the reg- 
ularity of three dimensional orbits of test particles. We know 
that and are piecewise functions with very different 
behavior when evaluated inside and outside the ring (fig- 
ures [2] and |3]) . The piece corresponding to the inner region 
describes the linear stability of the ring, whereas the one 
corresponding to the outer zone is associated to the regular- 
ity or chaoticity in three dimensional motion (NDCO). We 
showed that each one of these pieces tend to be joined at the 
outer edge by increasing the central body mass (figure [5]). 
This joining takes place in such a way that and i/^ tend 
to be positive valued functions throughout the equatorial 
plane. This means that more and more stable configurations 
lead to a gradual decrease of chaoticity in three dimensional 
orbits, and the phase-space structure associated with them 
tends to be of the type of regular toroids. 

It is worth to point out that the most relevant contri- 
bution to the stability in the configurations studied here, is 
provided by the mass of the central body. We can conclude 
that the increment of the monopolar term in the central 
body, rather than the decrease in the quadrupole moment, 
favors the existence of isolated islands in the phase space. 
The mass ratio and the size of the ring affects in the same 
way the stability of circular orbits (inside and outside the 
ring) and the regularity of three-dimensional motion (both 
DCO and NDCO). The fact that for DCO there exist a vari- 
ety of regular islands throughout the phase-space, is relevant 
to confirm that the circular orbits supporting the ring struc- 
ture are rather insensitive to small perturbations. Likewise, 
the above situation also holds for outer equatorial orbits and 
the associated three-dimensional motion (orbits without disc 
crossings) favoring the possible formation of new structures. 
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